Origin of negative differential thermal resistance in nonlinear systems 
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Negative differential resistance in electronic conduction has been extensively studied, but it is 
not the case for its thermal counterpart, namely, negative differential thermal resistance (NDTR). 
We present a classical Landauer formula in which the nonlinearity is incorporated by the self- 
consistent phonon theory in order to study the heat flux across a chain consisting of two weakly 
coupled lattices. Two typical nonlinear models of hard and soft on-site potentials are discussed 
respectively. It is shown that the nonlinearity has strong impacts on the occurring of NDTR. As 
a result, a transition from the absence to the presence of NDTR is observed. The origin of NDTR 
consists in the competition between the temperature difference, which acts as an external field, and 
the temperature dependent thermal boundary conductance. Finally, the onset of the transition is 
clearly illustrated for this model. Our analytical calculation agrees reasonably well with numerical 
simulations. 

PACS numbers: 05.70.Ln, 44.10.-l-i, 05.60.-k 



I. INTRODUCTION 

In the linear response regime, transport processes such 
as transport of mass, momentum or energy can be de- 
scribed by hnear laws of the form 



3 = t^F, 



(1) 



where j and F stand for the generahzed flux and force, 
respectively, and /i the transport coefficient. In other 
words, the flux j is a hnearly increasing function of the 
external field F, which is a well-known characteristic of 
Fick's law for mass transport. Ohm's law for electron 
transport and Fourier's law for heat transport. As the 
field F becomes too strong, the system may come into 
the nonlinear response regime, where the linear relation 
^ is no longer valid since the transport coefficient fi be- 
comes itself field dependent. As a result, an interesting 
phenomenon, i.e. negative differential resistance (NDR) 
may take place in a system in the strong-field regime 
where the flux counter-intuitively decreases as the exter- 
nal field increases. 

Since the pioneering observation in the tunnel diode 
by Esaki [Ij, NDR has been extensively studied for the 
electronic transport, which led to widespread practical 
applications in modern electronics [2]. It is still an ac- 
tive research topic to date, particularly at the atomic 
scale However, its counterpart in the heat conduc- 

tion problem, namely the negative differential thermal 
resistance (NDTR) has been much less studied. NDTR 
effect has been noticed in the studies of asymmetric heat 
conduction (see Fig. 1 in Refs. [5] and [6,]), where it has 
been shown to be critical to design a thermal diode with 
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an enormous rectification factor. It has also been shown 
that NDTR is crucial for a correct functioning of lat- 
tice models of thermal transistor 7] and thermal logic 
gates t^j. It is already known that NDTR effect can be 
qualitatively explained in terms of the overlappingof the 
vibrational spectra of the interfacial particles 0,9- On 
the other hand, is has been recently shown that in the 
model presented in 7], the NDTR effect will gradually 
disappear as the system size increases or the properties 
of the interfacial coupling change [l^ . For a clear un- 
derstanding of the mechanism underlying NDTR effect, 
it is thus imperative to comprehend from a quantitative 
point of view the necessary conditions for the occurring 
of NDTR effect. 

As far as we know, a rigorous analytical approach to 
study heat conduction in a non-integrable lattice system 
at the nonlinear response regime has been so far unavail- 
able. One usually has to rely on numerical simulations. 
In the present study, we will develop a phenomenological 
approach, in line with that in Ref. |6]] , to study heat flux 
through an "interface" between two weakly coupled an- 
harmonic segments. We study two typical models, which 
have hard and soft anharmonicity respectively. The the- 
oretical calculation based on the presented Landauer-like 
formula yields results in reasonable agreement with the 
numerical simulation. We will show that the intrinsic 
nonlinearity of the system is necessary for the occurring 
of NDTR. It is further iUustrated that NDTR does not 
always occur in the presence of nonlinearity. A transition 
from the absence to the presence of NDTR with the in- 
crease of the nonlinearity is illustrated for both models. 
A simple but physically appealing mechanism is proposed 
in order to explain the origin of NDTR in the nonlinear 
systems. Our study of NDTR provides indications of 
possible applications such as the construction of thermal 
devices. 
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II. THEORETICAL APPROACH 

We study in this work the stationary heat current 
across a chain consisting of two weakly coupled lattices, 

H^H+ + ^{x,^xof + H_, (2) 

where the Hamiltonian for the left and right segments 
are given by 

H+= ^ ^ + U^.+i-xO' + U+{x.) (3) 



and 



i=-JV/2+l 



N/2 2 , 

i=l 



(4) 



respectively. U±{x) represent the onsite potential that 
will be specified below. Two heat baths with tempera- 
tures T_|. and T_ are connected to the extremities of the 
left and right segment, respectively. Note that NDTR ef- 
fect have been so far investigated only in spatially asym- 
metric models However, we will show in this study 
that NDTR can also take place in a spatially symmetric 
model, i.e. for U+{x) = U^{x). 

In the case where the coupling Kint between the left 
and right segments is weak, the two segments will achieve 
two nearly equilibrium states at temperatures and 
T_, respectively. Their vibrational motion can thus be 
approximately described according to the self-consistent 
phonon theory (SCPT) [1, [Tl| - [l3j with effective harmonic 



Hamiltonians H^^'^ and H';^' of the form 



r(0) 



H 



(0) 



+ - Y + 2(^^+1 

i=-N/2+l 



\2 , J+ 2 



(5) 
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FIG. 1: Heat flux j as a function of T- for the symmetric 
harmonic chain (A = 0, /lo = /ho = /o). Here Kint = 0.05, 
r+ — 1. The system size is A'^ = 64 for the simulation. The 
linear dependence of j on T- implies that the transmission is 
temperature independent. 



for a given measurable ^(u), where /3± — {kBT±)~^. 
Note that according to Eq. ([51), f± is temperature de- 
pendent. The derivation of Eq. (O can be found in the 
appendix of Ref. [l^. The renormalized normal- mode 
frequencies of phonons in each segment can then be writ- 
ten as uj±{q) — y^4sin^ (9/2) + /±, where q stands for 
the wave vector. 

According to the Khalatnikov theory [14], the heat 
flux can be regarded as the propagation of plane waves 
(phonons) with various frequencies. Within this ap- 
proach, the heat flux through the system composed of 
the segments ^ and ^ can be determined by [Tsj 



J 



kB{T+-T^) 

271 



T{ijj) dw. 



(9) 



2 1 f 
tt(0) ST Pi It \2 , /- 



(6) 



with Ui = Xi — (xi) = Xi — rj. The effective force con- 
stants f± are determined by numerically solving the self- 
consistent equations 



/± = 2 



d{x-^)f 



(7) 



where Wmm and ujmax correspond to the boundaries 
of the overlap band of left and right phonon spec- 
tra, that is uJrnin = Hiax { y^Tf , \/J^} and 
min {y^ A + ■y/4 -|- /_ } . T{uj) is the phonon transmis- 
sion probability through the interface. It is worth men- 
tioning that Eq. (|9]) is similar in form to the celebrated 
Landauer formula 



1 

2^ 



T{U) [77+ {lJ,T+)-7J^{uJ, r_ )]nu duj, (10) 



Here denotes the thermal average with respect to ^liere r,± = {6^^ ^-1) are the Bose-Einstein distri 

the trial harmonic Hamiltonian at the correspond- 
ing temperature T± and it is defined by 

/ A{u) exp (-/3±i7i°^ (u)) du 



(A(u))L°) 



/exp(-/3±i7f (u))du 



(8) 



bution functions. The Landauer formula ([TU)) . originated 
from the study of electron transport [l^, describes the 
ballistic transport of phonons in quantum systems |l7{ . 
Considering the high temperature limit (classical limit) 
where r]± = {(3±hu;)~^, Eq. reduces to Eq. ©. Note 
that the quantum constant h is cancelled automatically 
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FIG. 2: Heat flux j as a function of T- for the 0* model. Left: 
A = 0.8; Right: A = 5. In both cases, Kim = 0.05, T+ = 1 
and TV = 64 for the numerical simulation. 



in the classical limit. Thus our equation (jH]) can be con- 
sidered as the classical form of the traditional Landauer 
formula. 

To find the transmission coefficient, we consider a 
plane wave incident from the left, which is partly refiected 
by the interface with amplitude R and partly transmits 
across the interface with amplitude T into the right seg- 
ment [6] . The displacement of the ith particle from the 
equilibrium position can be written as 



{e^''' + Re 

rpglq{i-l)-lu]-t 



Iki\^ — Iiu^t 



if i < 

if i > 1 



(11) 



where / is the imaginary unit. Thus the motion of the 
interface particles can be described by the following equa- 
tions 

~uj\uq = u_i + K.„itui - (1 + Kint + f+)uQ, (12a) 



-W^Ul = U2 + KintUo - (1 + Kint + f-)ui 



(12b) 



If an acoustic matching condition lu^ = uj^ = cj is 
satisfied, the solution of Eq. gives the transmission 
probability of Eq. (j9l) in the form 



r(c.,/+(T+),/_(T_) ) = l-|i?2 
Ci(l-2if„t) + C3ifL' 



(13) 



where 



Csii^) = (Ci + C2)/2 + 2^2 _ _ /_ . 



(14a) 
(14b) 
(14c) 

Thermal transport is inhibited, i.e. T = if the phonon 
bands of the two segments are mismatched. 

In what follows, we will illustrate the results of the 
analytical calculation based on Eq. As a compari- 
son, non-equilibrium molecular dynamics (NEMD) sim- 
ulations are performed by applying Langevin heat baths 
at the two extremities of the chain. The heat fiux, whose 
definition can be found in Ref. [isj . is averaged over a 
long enough time so that the system reaches the steady 
state regime, at which the local heat flux is constant along 
the chain. During the simulations, T_|_ was fixed and the 
temperature difference AT = T+ — T_ was changed by 
changing T_ . 




FIG. 3: (f* model: The turning point T*, at which the heat 
current exhibits a maximum, as a function of A. Nonzero 
T* indicates the presence of NDTR. The transition from the 
absence to the presence of NDTR occurs at Ac f« 1. 



A. (j) model 

Let us first consider the cj)'^ model, which is a typical 
bounding potential of "hard" anharmonicity , 



TT f \ fo 2 : 4 
U±{x) = —X^ + -X^ 



(15) 



The effective force constants /± are determined by nu- 
merically solving the self-consistent equations 



/± = /o 



SAfc^Ti 



4/i 



(16) 



Note that j± regularly increases with increasing temper- 
ature. 

Before discussing the nonlinearity effect, we will ap- 
ply the above analytical analysis to the harmonic model 
A = 0. In this case, the transmission (|13p is tempera- 
ture independent since /± = /q. Fig. [T] illustrates the 
heat flux j as a function of r_ for the harmonic chain 
for several values of the harmonic constant /q. By in- 
specting the figure, we first notice that in the ballistic 
case, the simulation agrees well with the analytical re- 
sult that follows from the classical Landauer formula (|9]) 
(ks = 1). Furthermore, it is seen that j increases lin- 
early when r_ decreases, that is when the temperature 
difference increases. As is expected, NDTR cannot be ob- 
served in the harmonic model since there does not exist 
any nonlinear response mechanism. 

Furthermore, Fig. [5] shows that if the nonlinearity is 
present inside each segment (A 7^ 0), the formula ^ still 
works reasonably well. By comparing the left and right 
plots of Fig. [21 it is seen that a transition from the absence 
to the presence of NDTR occurs as the nonlinearity A 
increases. We also present in Fig. [3] the evolution of the 
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turning point T* at which NDTR effect manifests itself 
as a function of the nonUnearity. The on-set of NDTR at 
Ac ~ 1 is clearly shown. 

Note that Eqs. (fTB|) and ^ implies the following scal- 
ing relation 



(17) 



where s is an arbitrary scaling constant. The same scal- 
ing property can be obtained from the equations of mo- 
tion of the model ©(see [19]). Eq. (H?]) indicates the 
equivalence between the temperature and the nonlinear- 
ity. Thus a similar transition for fixed A can be observed 
as T+ increases, which will be verified both from numer- 



ical simulations and our analytical approach. 
NDTR takes place if A > Ac (or T+ > Tc). 



B. On-site Morse model 



In fact, 
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FIG. 4: Effective force constant / as a function of the tem- 
perature for the model (jlSp . Inset: The critical temperature 
Tc, above which / vanishes, against the nonlinearity D. 



Now we consider a model which consists of two weakly 
coupled symmetric nonlinear lattices with an on-site 
Morse potential given by 

U±{x) ^ D[eii-p {~ax) ~ l]"^ . (18) 

Model (|18p was introduced in order to investigate the 
DNA dcnaturation process [13] • The anharmonicity of 
the model is "soft" since the Morse potential is bounded 
for X — >■ oo. 

The effective force constants f± are obtained by nu- 
merically solving the following self-consistent equations 

f^ = 2a^Do.J^^^M^]. (19) 

It should be noted that f± decrease as the temperature 
increases for a soft anharmonic potential like Eq. ([TS]). 
as shown in Fig. 21 This is different from that of a 
model with a hard anharmonicity like Eq. (1151) , for which 
f± monotonically increases as the temperature increases. 
One can see from Fig. |3]that there exists a critical tem- 
perature Tc, above which the force constant / vanishes. 
It means that the on-site potential can be neglected once 
the thermal energy of the particles overcomes the poten- 
tial energy and then the system behaves like a harmonic 
chain. The inset of Fig. |4] shows that the critical tem- 
perature increases with the nonlinearity of the system, 
which reflects the fact that the deeper the potential well 
is, the larger is the thermal energy needed to overcome 
the potentiel barrier. 

We then use Eq. © to compute the heat fiux and com- 
pare the analytical result with numerical simulations, as 
shown in Fig. [SJ One can see that since the nonlinearity 
D is weak, the heat flux increases monotonically with in- 
creasing AT . Nevertheless, NDTR occurs as D becomes 
large enough. Like cj)'^ model. Fig. [S] indicates that as 
the nonlinearity increases, the soft potential model also 
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FIG. 5: Heat flux j as a function of T- for the Morse model. 
Left: D = 5 X 10"^; Right: D = 15. In both cases, a = 
0.426, Kint = 0.05, T+ = 20 and iV = 64 for the numerical 
simulation. 



exhibits a transition from positive differential thermal re- 
sistance (PDTR) to NDTR. 

A simple scaling relation like Eq. (flTl) for the nonlin- 
earity and the temperature is inexistent for the Morse 
model. However, increasing the nonlinearity is qualita- 
tively equivalent to increasing the temperature. Thus a 
similar transition from PDTR to NDTR as the temper- 
ature increases can be expected. In Fig. [51 we calculate 
the scaled turning point r* = T* /T+, at which the heat 
flux exhibits a maximum. Non-zero r* indicates the pres- 
ence of NDTR. The analytical result shows the on-set of 
NDTR at r+ sa 11 for the dissociation threshold fixed at 
D = 15. 

However, one should note that the SCPT fails in the 
vicinity of the melting transition for the present soft po- 
tential model, as pointed out in Ref. The reason 
for the failure of the variational approach is the fact 
that the half-bounded potential (fTS]) cannot be simply 
approximated by the trial harmonic potential with a 
temperature-dependent force constant. This pecularity 
prevents us from modelling the crucial role of the non- 
linearity in this regime using the SCPT. This point will 
be discussed in detail in the next section. Although it is 
clear that the quantitative agreement with the simulation 
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FIG. 6: Morse model: The scaled turning temperature 
T* /T+, at which the heat current exhibits a maximum, as 
a function of T+. Here D = 15,a = 0.426. 

result is poor, we emphasize that the present approach 
goes far beyond the traditional Landauer approach in its 
ability to characterize the nonlinear response regime, for 
which a transition from PDTR to NDTR is illustrated at 
least in a qualitative manner. 

III. PHYSICAL MECHANISM 

The results presented so far give rise to the follow- 
ing question: what is the origin of NDTR in the above 
models? We will show that the classical Landauer equa- 
tion © yields a simple and intuitive explanation. One 
should note that the temperature discontinuity at the 
virtual interface(the site xq) indicates the existence of 
the thermal boundary resistance (or conductance, see 
Ref. [21'!), and it plays a crucial role for the heat conduc- 
tion in our weakly-coupled model. Defining the effective 
thermal boundary conductance by 

£'^'W{u)duj, (20) 

Eq. ^ can be rewritten as 

j = ATa, (21) 

which is similar in form to Ohm's law for electrical con- 
duction. The simple relation (pij) suggests that there 
exists mainly two contributions to the temperature de- 
pendence of the heat current for a two-segment system. 
The first contribution comes from the temperature dif- 
ference AT, which acts as an external thermal force and 
yields the regularly increasing behavior of j with decreas- 
ing T_. The second contribution is due to the thermal 
boundary conductance a. One can see from Fig. [7] that 
unlike AT, a is an increasing function of T_. The widen- 
ing of the overlap band of the vibrational spectrum of 
segments L and R, or Aw = uJmax — i-^min, is mainly 
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FIG. 7: (j)'^ model: The thermal boundary conductance a (left 
plot) and the heat flux j (right plot) against T- for different 
values of A. Lines 1, 2, 3 correspond to A = 1, 4, 9, respec- 
tively. One can see that nonlinearity A has an important effect 
on the behavior of the heat flux due to the interchange of the 
dominant role in Eq. (pTj) . Here T+ = 1, Kim = 0.05. 



responsible for this increasing behavior of the thermal 
conductance. Consequently, the origin of NDTR effect is 
basically the competition between the growing "external 
field" AT and the diminishing overlap band Aw as T„ 
decreases. NDTR thus occurs below the temperature T* 
at which the opposite behavior of both contributions ex- 
actly compensate each other and it takes place if and only 
if a becomes dominant for T_ < T+ . Fig. [7] shows the 
apparition of NDTR effect as one considers different non- 
linearity parameters A = 1, 4, 9 for the 0"* model. One can 
note that for segments with A large enough, j vanishes as 
T_ decreases due to the mismatch of the phonon bands. 
We also plot in Fig.|S]the thermal boundary conductance 
a and the corresponding heat flux j for the Morse model, 
which displays a similar behavior. We can thus conclude 
that the proposed mechanism for the occurring of NDTR 
is valid for both hard and soft models. For the harmonic 
system, a is exactly temperature independent, leading to 
the linear behavior observed in Fig. [1] 

One should note that the curve 2 of Fig. [5] exhibits 
a jump at T_ = Tc . For T_ > Tc, the effective force 
constant / vanishes and the phonon frequency becomes 
temperature independent. The system thus behaves like 
a harmonic chain, characterized by a temperature inde- 
pendent thermal boundary conductance a and a linear 
behavior of the heat flux. The occurrence of the jump 
and the exact linear behavior of the heat flux as T_ > Tc 
are inconsistent with the numerical simulation. As dis- 
cussed in the last section, this artificial result lies in the 
incapability of the SCPT to model the transition at the 
vicinity of Tc as shown in Fig. 2) Since Tc increases with 
D, the artificial jump disappears provided T-|_ < Tc{D), 
which is illustrated in the right plot of Fig. [S] 



IV. SUMMARY 

In summary, we presented a classical Landauer formula 
to study NDTR effect in typical lattice models. It was 
shown that NDTR cannot occur in a harmonic lattice, for 
which the linear relation ^ is generally valid no matter 
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FIG. 8: Morse model: The thermal boundary conductance a 
(left) and the heat flux j (right) as a function of T_ for differ- 
ent values of D. Lines 1, 2, 3 correspond to D = 0.001, 5, 18, 
respectively. Here T+ = 20, Kint = 0.05. 

how large the temperature difference is. In the presence 
of anharmonicity, one can observe a transition from the 
absence to the presence of NDTR as the nonlinearity is 
increased for both hard and soft potentials. The NDTR 
effect is basicahy due to the competition between the in- 
creasing behavior of the external field and the decreasing 
behavior of the effective thermal boundary conductance 
of the chain. The transition in Figs. [3] and [6] indicates 
that NDTR may be controlled by adjusting the param- 
eters of the system or the temperature of heat baths, 
whose utility for nanoscale applications is clear. 

It is imperative to clarify why we are allowed to apply 
such a seemingly ballistic transport formula (O to cal- 
culate the heat flux through a nonlinear system in the 



strong field regime. For the particular model presented, 
even though the whole system is at strong external field, 
each segment is still close to its corresponding equilib- 
rium state thanks to the weak coupling, and can thus be 
approximately described by SCPT. It should be stressed 
that the analytical estimation, which is based on the lo- 
cal thermal equilibrium of the segments, holds only if the 
coupling Kint is weak enough. For strong coupling, both 
segments get far from thermal equilibrium and SCPT can 
not be applied anymore to deal with the nonlinearity. 

Finally, it is worth giving a comment about the rela- 
tion between asymmetric heat conduction and NDTR. 
Note the following two facts: 1) Asymmetric heat con- 
duction results from the intrinsic spatial asymmetry of 
the model, which is not necessary for the occurring of 
NDTR as shown in this study; 2) One can observe asym- 
metric heat conduction without the occurring of NDTR 
as long as the applied temperature difference is moder- 
ate. It seems that the NDTR, as a field-induced effect, is 
neither a sufficient nor a necessary condition for thermal 
rectification. 
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